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Abstract 


In this paper, we are interested in a general equation that has finite speed of propagation 
compatible with Einstein’s theory of special relativity. This equation without external force 
fields has been derived recently by means of optimal transportation theory. We first provide 
an argument to incorporate the external force fields. Then we are concerned with comparison 
and maximum principles for this equation. We consider both stationary and evolutionary 
problems. We show that the former satisfies a comparison principle and a strong maximum 
principle. While the latter fulfils weaker ones. The key technique is a transformation that 
matches well with the gradient flow structure of the equation. 


1 Introduction 

1.1 Motivation 

This paper is concerned with comparison and maximum principles for the following equation 


dtu = div[uV(/J*(VU)] + div[tiV(/J*(V logu)], 


( 1 ) 


In this equation, the spatial domain is IR"^; the unknown is tt : [0, T] —div denotes the diver¬ 
gence operator; the function V : R is given, and finally (p* is the Legendre transformation of 

a given function ip G A typical example of Eq. ([T]) is when p is the so-called relativistic 

cost function 
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where c is a given constant which has the meaning of the speed of light. In this case, 

X 


^*c{x)=c + , and Vip*{x) = 


1 + 


|x|2 


and Eq. © becomes the well-known relativistic heat equation (with external force fields) 


dtu = div 


uVV 


1 + ^ 


div 


uVu 




In general, in this paper, the function ip will be of general form 

c(|a;|), if |a;| < c. 


ip{x) = 


- 1-00 if Iccl > c. 


where c is a given function on the real line. 


(3) 


(4) 


Let us first review relevant literature about the origin of Eq. ©• It is well-known that the 
classical heat equation, dtU = Au, satisfies a strong maximum principle which implies infinite 
speed of propagation. However this is incompatible with Einstein’s theory of special relativity 
that the speed of light is the highest admissible free velocity in a medium. Eq. which belongs 
to a larger class of the so-called flux-limited diffusions, have recently been introduced to fulfil 
Einstein’s theory. In particular, the relativistic heat equation (i.e., Eq. ([3]) without external field 
E = 0) was proposed first by Rosenau [23j and was derived later by Brenier by means of optimal 
transportation theory. He formally showed that the relativistic heat equation is a gradient flow of 
the Boltzmann entropy with respect to a Kantorovich optimal transport functional associated to 
the relativistic cost (I2|). In [^, the authors proved this rigorously and generalized to the general 
cost in (ID) obtaining Eq. © with V = 0. Eq. ([3]) with a non-zero force field has also appeared 
in other contexts. For instance, to be able to describe the velocity profile of charge carriers along 
submicrometer electronic devices under high fields, in [27], the author obtained Eq. dSj) using the 
method of continued-fraction expansion. As another example, in the context of the flux-limited 
chemotaxis model, in [^ the author derived a similar equation as (I3|) by optimizing the the flux 
density of particles along the trajectory induced by the chemoattractant. Eq. © with a non-zero 
external force field is a generalisation of Eq. ([Sj) to the general cost function (|3|). 

Eq. (|T|), and in general flux-limited diffusion equations, play an important role both in 
theoretical research and in practical applications. From a practical point of view, these models have 
been used widely, for instance, in plasma physics Hi, in the modelling of the charge transport in 
microelectronic devices [IHIIIH], and of multicellular growing systems mm- Besides their practical 
usage, they also create interesting mathematical problems due to the nonlinearity and discontinuity 
of the function ip. To solve these problems, many authors have established intriguing connections 
between different fields of mathematics such as nonlinear partial differential equations, optimal 
transportation and differential geometry, see e.g., EiisiiiniiH] and references therein. 
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1.2 Aim of the paper 

The aim of this paper is to study comparison and maximum principles for Eq. If}. Com¬ 
parison and maximum principles have been the subjects of intensive research, see the mono¬ 
graphs [MUZUll] and references therein for more information. These principles provide insightful 
information about partial differential equations such as propagation of speed and uniqueness. As 
mentioned in the previous section, Eq. © has finite speed of propagation; hence it certainly does 
not satisfy a strong maximum principle. Therefore, it is natural to ask whether it does fulfil a 
weak maximum principle. This is the focus of the present paper. We consider both the stationary 
and evolutionary problems. We show that the former satisfies a comparison principle and a strong 
maximum principle. While the latter fulfils weaker ones. It should be mentioned that in a recent 
paper [3T] , the authors address this issue affirmatively for the case of the relativistic heat equation 
(i.e. Eq. dS]) with E = 0). The present paper extends their results to Eq. (U]) with non-zero 
external force helds. As presented above, Eq. ([1]) covers a much larger class of equations that are 
of importance in both theoretical research and practical applications. 

We now summarize our assumptions and main results. 

1.3 Assumptions and main results 

Assumption 1.1. Throughout the paper, we make the following assumptions. 

(HI) V G C'2(1R'^), 

(H2) (f* G C'^(1R'^) is strictly convex (i.e., all eigenvalue of the Hessian are positive). 

We now describe our main results. Let Lu denote the operator on the right hand side of Eq. 
O, i.e., 

Lu = div[uV(/J*(VH)] -I-div[uV(^*(V logrt)]. (5) 

We consider both the stationary problem 

Lu = 0, (6) 

and the evolutionary one 

dtu = Lu. (7) 

The first theorem is a tangency principle m for the stationary problem. 

Theorem 1.1 (Tangency principle for the stationary problem). 

Suppose that H C IR'* is open and connected. Assume that u,v € (7^(12) and u,v > 0. If Lu > Lv 
and u < V in n and u = v at some point in LI, then u = v in LI. 


4 


We obtain the following strong maximum principle as a direct consequence of the tangency 
principle. 

Theorem 1.2 (Strong maximum principle for the stationary problem). 

Suppose that il C R‘^ is open and connected. Assume u G H (7(fl) and u > 0. 

(1) Additionally, assume that div(Vv3*(Vy)) < 0. If Lu > 0 in VL and u attains its maximum 
over n at an interior point, then u is constant within 0. 

(2) Similarly, assume in addition that div(V(/5*(VT)) > 0. If Lu < 0 in U and u attains its 

minimum over O at an interior point, then u is constant within LI. 

A comparison principle also holds true for the stationary problem. 

Theorem 1.3 (Comparison principle for the stationary problem). 

Suppose u,v G C'^(Ll) fl C(ri) with u,v > 0 and satisfy Lu > Lv in LI. If u <v on dLl, then u < v 
in LI. 

Next we deal with the evolutionary case. Let 12 C R'^ be an open set and T > 0 be given. We 
denote the space-time domain by LIt := LI x [0,r]. The boundary of LIt is = LIt \ LIt- Finally, 

C^d(57y) denotes the set of functions of {x,f) G LIt which is in space, and in time. 

Theorem 1.4 (Comparison principle for the evolutionary problem). 

Suppose that u,v G C'^'^{L1 t) H CILIt),u,v > 0 satisfying dtu — Lu < dtv — Lv in LI. Then 
max('u — v) = max(it — v). In particular, if u <v on Ft, then u <v in LIt. 

Qt Tt 

As a consequence of the comparison principle we obtain the weak maximum principle for 
evolutionary problem. 

Theorem 1.5 (Weak maximum principle for the evolutionary problem). 

Suppose that u G C^'^(LIt) H C(LIt) and u > 0. The following statement hold. 

1. Additionally, assume that div(V(/3*(Vy)) < 0. If dtU — Lu < 0, then a maximum of u in LIt 
is attained at the boundary Ft. 

2. Similarly, assume in addition that div(V(p*(Vy)) > 0. Ifdtu — Lu > 0, then a minimum of 
u in LIt is attained at the boundary Ft. 

The underlying idea of the proofs of the these theorems is to transform the operator L into 
a strictly elliptic operator Q that is independent of the function itself. The assertions for Q are 
deduced from well-known results on maximum/comparison principle for quasi-nonlinear partial 
differential equations in [m iizi US]- Then we obtain the corresponding results for the original 
operator L using the monotonicity of the transformation. 
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1.4 Organisation of the paper 

The rest of the paper is organized as follows. In Section[2J we summarize the derivation of the 
main equation. In Section [3l we prove the main theorems for the stationary problem. Section |4] is 
devoted to the evolutionary case. Finally, we conclude the paper and discuss future perspectives 
in Section [5] 


2 On the derivation of the main equation 


In this section, we derive Eq. m and prove that, under an appropriate assumption on the 
growth of V, it has an unique equilibrium solution given by the Gibbs measure. 

In a seminal paper |19| . Jordan, Kinderleher and Otto proved a remarkable result that the 
classical heat equation 

dtu = Alt (8) 

is a gradient flow of the Boltzmann entropy S{u) = f (u log u — u) with respect to the Wasserstein 
metric. This statement means that the heat equation can be obtained by computing the time 
discrete solution at the time step k as follows (now known as JKO-scheme) 

= argmin^g.p(Kd)S'(u) + u), (9) 

where the Wasserstein distance W 2 ^h is defined via 

Wlh{u^~\u)= inf / l{dx,dy), 

and passing to the limit as the time step h goes to 0. In the expressions above, P(IR^) is the set of 
all probability measures on IR'^ and r(^, denotes the set of probability measures on whose 
marginals are /i and v. To obtain the linear drift-diffusion equation 


dtu = div(uVF) J- Au, 


( 10 ) 


one simply replaces the Boltzmann entropy S by the free energy F = S + J Vu. Since then 
the paper m has embarked a lot of research that links many branches of mathematics together 
such as partial differential equation, optimal transportation, differential geometry and probability 
theory, see e.g., [nmnniiii], the monographs mus] and references therein for more information. 

In [3], Brenier formally derived the relativistic heat equation (i.e., Eq. ([3]) with V = 0), from 
the scheme (HI) replacing IE|^ by 








Wc = inf 
7er(u'“-i 


( 11 ) 
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where ipc is the relativistic cost function (ED. Later on, McCann and Fuel m proved this rigorously 
and extended to the general case with ip given in 0, obtaining the following equation 

= div[uV<p*(Vlogu)], (12) 

which is Eq. o where the external force is absent. 

In contrast to the linear drift-diffusion equation dlUD, it is not obvious how to incorporate 
the external force fields into Eq. (US- One possibility, as in the classical setting, is to replace 
the Boltzmann entropy S = J(ulogu — u) by the free energy F = S + f Vu. This amounts to 
substitute = log u in (IT^ by = log u + VE leading to the following equation 

dtu = div[MV(^*(logM -I- VE)]. 


For the relativistic cost, since Vv3*(z) = . ^ ^ the above equation becomes 


dtU = div 


u 


log u -I- VE 



I logii+VVP 

c2 


(13) 


Note that due to the nonlinearity of '\/ip*{z), this equation is different from Eq. (ED. In comparison 
with the classical equation dloD, the drift and the diffusion terms are separated since V(p*(z) = z 
is linear. We do not know whether the equation above has any physically/biologically meaning. 
Therefore, we take different approach. Before that we recalling the derivation of (ED. To obtain 
this equation, different methods have been exploited before in the literature. For instance, m 
used continued-fraction technique to derive (ED from the hierarchy of hydrodynamic equations. 
As another example, in the context of the flux-limited chemotaxis modelling, in [B] the authors 
obtained it by optimizing the the flux density of particles along the trajectory induced by the 
chemoattractant. Note that in Eq. 0, the drift term is equal to div[MVv3*(VE)]. Motivated by 
this observation, we now include the external force into (1121) simply by generalising the above term 
to div[uV(/7*(VE)] and adding it to (fT^ . This results in the following equation 


dtU = div[uV(/j*(VE)] -I-div[itV(/j*(V logu)], (14) 

which is exactly Eq. o that we started with. When ip = ipc, we do get Eq. ED- One important 
property of both Eq. m and the classical linearly drift-diffusion equation is that the Gibbs 
measure Z~^ exp(—E) is the unique equilibrium solution. In Lemma l2.II below, we prove that this 
property is preserved in Eq. O- As an additional observation, since iPci^) —z as c —> oo, at 
least formally we recover the classical drift-diffusion equation from Eq. m as c —>• oo. It would 
be interesting to investigate whether one could derive Eq. (ED solely from the optimal transport 
theory as in |20) and recover its classical counterpart as c goes to infinity . We provide more 
discussion about this in Section [S] 

We now show the following. 
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Lemma 2.1. Eq. (fTOl) has an equilibrium solution which is given by 

Ueq{x) = CexY>{-V{x)), (15) 

for some constant C. If V satisfies that f exp(—V(x)) dx < oo, Ueq can he normalized to become 
the Gibbs (probability) measure. 

Proof. The proof will use the following properties of 

(PI) is odd, i.e., V(p*{z) = —V(p*{—z). This is because p*{z) = f{\z\) for some function /. 

(P2) \7ip* is injective, i.e., if S/ip*{zi) = V(p*(z 2 ) then = Z 2 . This is because (p* is a convex 
function that implies that is monotone. 

We recall that u is called an equilibrium solution of Eq. © if it annihilates the total flux. Now, 
suppose that Ueq is an equilibrium solution, i.e., 

Vip*(yv) + Vip*{VlogUeq) = 0. 

Hence by (PI) 

yip*iyiogueq) = -V(^*(VP) = V(p*(-vp), 

and then by (P2), 

V log Me, = -VP, 

which implies that Me, = C'exp(—P) for some constant C. If f exp(—V(x)) dx < oo, there exists 
only one value of C, namely C = (f exp(—P(a;)) dx) such that Me, is a probability measure on 
R'^. This is exactly the Gibbs measure. 

However, we note that Eq. o is still meaningful without this extra condition. This is also 
true for the classical drift-diffusion equation. In addition, its Wasserstein gradient flow structure 
is rigorously well-defined without the extra assumption on the growth of P, see [19j . □ 


3 The stationary problem 

In this section we deal with the stationary problem O, i.e., 

Lu = 0, (16) 

where the operator L is defined in ([5]). 

We consider a change of variable w = logM, and transform Eq. dUD into the following one 


Qw = 0 , 


( 17 ) 


where the operator Q is defined by 


Qw = div[V(/j*(Vu;)] + Vw • V>p*{Vw) + Vw ■ Vip*{VV) + div Vv3*(Vy). (18) 

The transformation plays a crucial role in the present paper. We provide more discussion about 
it in Remark 13.51 

We now show that the stationary problem (EZD satisfies a tangency principle, a strong maxi¬ 
mum principle and a comparison principle. Let us begin with the tangency principle. 

Proposition 3.1 (Tangency principle for the operator Q). Suppose that C R'^ is open and 
connected and that w,w' S If Qw > Qw', w < w' in O, and w = w' at some point in fl, 

then w = w' in Q. 

Proof. We first show that Q is a strictly elliptic operator. Indeed, since 
div[V^-(Vu,)| = = VV(V») : 

the operator Q is equal to 

Qw = VV*(Vw) : -k Vw ■ Vv3*(Vw) -k Vw • V</j*(VV) -k div V(/5*(VV). (19) 

We can write Q in the form of a quasilinear operator studied in the by-now classical monograph 
[TTI Chapter 10] 

d 

Qw = ^ a*'^ (a;, w, Vw)9j^j w-|-6(a;, w, Vw), (20) 

*.i=i 

where the coefficients of Q, namely the functions (x, z,p), i,j = l,...,d and b{x, z,p) are defined 
on fl X R X Note that the coefficients are independent of z: b depends on x,p, while a depends 
only on p. More precisely, 

aix,z,p) = y^ip*ip), b{x,z,p)=p-yip*{p)+p-Vip*iyv) + dwVip*iyv). ( 21 ) 

Since (p* is strictly convex, the matrix 'V'^(p*{z) is positive definite. Therefore, Q is a strictly 
elliptic operator. Furthermore, since the coefficients of Q are independent of w and continuous 
differentiable in Vw and V^w, the assertion of the proposition is obtained by applying the tangency 
principle for nonlinear elliptic operators in [221 Theorem 2.1.3]. 

□ 

The tangency principle is a strong statement. It allows us to obtain the following strong 
maximum/minimum principle for the operator Q. 

Proposition 3.2 (Strong maximum principle for the operator Q). 

Suppose that C R'^ is open and connected. 
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(1) Assume that w G fl C{^1) and that div(V(^*(Vy)) < 0. If Qw > 0 in fl and w attains 

its maximum over at an interior point, then w is constant within H. 

(2) Similarly, assume that w' G nC'(ri) and that div(V(/3*(Vy)) > 0. If Qw' < 0 in Gl and 

w' attains its minimum over Gl at an interior point, then w' is constant within Gl. 

Proof. (1) In ProDOsition l3.ll we take w' = M, which is the maximum attained by w. We have 

Qw' = div(V(/5*(VP)) <0<Qw. 

In addition, w < M, and u = M is attained at some point in Gl. The strong maximum 
principle is then followed straightforwardly from the tangency principle. 

(2) Similarly the minimum principle is proved by taking w = m', which is the minimum attained 
by w'. 

□ 

Remark 3.3. In the case of the relativistic heat equation Q, with (for simplicity we put c = 1) 

<P:(x) = ^/T+W - h ( 22 ) 


we have 


\/l + |a;p 


1 



XiXj \ 

i + \x\y' 


Hence, we obtain 


Qw 


V^ip*{Vw) : V^w + Vw ■ Vip*{Vw) + Vw ■ V(p* (VP) + div S/cp* (VP). 
A'u;(l + |Vwp) — V^rc : (Vrc (8> Vrc) |Vwp Vw • VP 

(1 + |2|2)i ^ + ^ '^/TTWW^ 

AP(1 + |VP|2) - V2p : (VP ® VP) 

(l + |VP|2)f ■ 

d 

a''^{x,w,Vw)dfjW + b{x,w,Vw), 

i j'=i 


where 


a^j{x,z,p) 
b{x,z,p) = 


= (vV*(p))., = 

\P? , 


V^TW 

p- VP 


^TTW Vi + ivpp 


+ 


. __PGP3_ 

1 + bP 

AP(1 + I VPp) - V2p : (VP (g VP) 
(1 + |VP|2)i 


It is interesting to note that in this case, the principle part of Q is exactly the mean-curvature 
operator. This raises some interesting questions, which we comment further in Section [S] 
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In particular, in the one-dimensional setting d = 1, we have 


div(V</5*(Vy)) 


V"{1 + V'^) - V"V'^ 
{1 + V'^)i 


V" 

(l-hF'2)f 


(23) 


Then the condition on V in the maximum/minimum principle in Proposition [X2] is easily verified: 
div V(^*(Vy) < 0 (> 0) if and only if V is concave (convex, respectively). □ 

We proceed with a comparison principle for the operator Q in the following proposition. 


Proposition 3.4 (Comparison principle for the operator Q). Suppose w,w' G fl C(fl) and 

satisfy Qw > Qw' in ft. If w < w' on dU, then w < w' in Q. 


Proof. Similarly as in the proof of Proposition 13.11 this proposition is obtained by applying the 
comparison for nonlinear elliptic operators in [53 Theorem 2.1.4]. □ 


We are now ready to provide the proofs of the main results for the original stationary problem. 

Proofs of Theorem \1.1\. Theorem IJ.2I and Theorem I J. ,91 


Since the exponential transformation u = exp(ri;) is positive and monotone, the statements of 
Theorem o Theorem O and Theorem 11.31 respectively follow from those of Proposition 13.11 
Proposition [33] and Proposition [331 □ 

We now discuss more on the transformation between u and w. 

Remark 3.5. As one can recognise from the proofs in this section, the transformation w = logit 
has two advantages. Firstly, it makes the operator Q depends only on Vic and not on w itself. 
This enables us to obtain maximum and comparison principles for the operator Q by applying 
known results for quasilinear elliptic operators. Secondly, it is positive and monotone, therefore, 
we achieve results for P straightforwardly from the corresponding ones for Q. 

The transformation is clearly hinted by the form of Eq. O- Tracing back from its derivation in 
Section|3 it is the gradient flow structure of the Boltzmann entropy with respect to the Kantorovich 
optimal transport functional associated to the general relativistic cost that gives rise to that form. 
In other words, the transformation matches perfectly well with the structure of the equation. 

The exponential transformation has been frequently used in the literature for the study of 
nonlinear partial differential equations (PDEs). However, note that the exponential transformation 
in this paper, is different from the well-known Cole-Hopf transformation. The latter is used to 
convert a nonlinear into a linear PDE. While we still get a nonlinear PDE in the former. If ip is 
the classical cost function, ip{x) = the transformation in this paper is actually inverse of the 
Cole-Hopf one. □ 
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4 The evolutionary problem 

In this section, we prove the main results, Theorem II.41 and Theorem II.51 for the evolutionary 
problem [7] 

dtu = Lu. 

Using the transformation as in the stationary case, the above equation is transformed to 

dtw = Qw. (24) 

The strategy now is analogous to the stationary case in the previous section. We prove comparison 
and maximum principles for the evolutionary problem associated to the operator Q using known 
theories in the literature, and then obtain the corresponding results for the original problem for 
the operator L using the positivity and monotonicity of the transformation. 

We recall that the definition on the spacetime domain and functions are given in the In¬ 
troduction. We first consider a comparison principle for the evolutionary problem associated to 

Q. 

Proposition 4.1 (Comparison principle for the evolutionary problem for Q). 

Suppose that w, w' G H C{flT) satisfying dtW — Qw < dtw[ — Qw' in fl. Then max(w — 

w') = max(r(; — w'). In particular, if w < w' on Ft, then w <w' in VLt- 
Tt 

Proof. The proof follows the procedure as in EH Theorem 3.3], which is in turn mainly adapted 
from [HI Theorem 10.1]. We rearrange the inequality in the assumption to obtain 

dt{w — w') — {Qw — Qw') < 0. (25) 

Set b{p) = p ■ V(p*{p) +p ■ V(/3*(VU). We have 

Qw — Qw' = a(Vr(;)V^w -I- b{yw) — a{Vw')'\/^w' — b{Vw') 

= a{Vw){V^w - V^w') + {aiVw) - a{Vw'))y‘^w' + b{Vw) - b{yw'). 

We set 2 = w — w', and define the linear operator Qz = a{x, t)V^z + b{x, t)Vz, with 

a{x, t) = a(Vic), 

b{x, t)yz = {a{Vw) — a(Vw'))V^w' + b(Vw) — b(Vw'). 

Note that the existence of the function b is due to the mean value theorem. Then we obtain 
Qw — Qw' = Qz, and (ITS!) becomes 

dtz - Qz <0, 

in rix- In addition, z < 0 on Tt. By the parabolic weak maximum principle m Theorem 8, 
Section 7], we have that z < 0 in fix, i.e., w <w' in VLt- This completes the proof. □ 
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As a consequence of the comparison principle we obtain the weak maximum principle for 
evolutionary problem. 

Proposition 4.2 (Weak maximum principle for the evolutionary problem for Q). 

The following statements hold. 

(1) Suppose that w € (fix) (fir) and div(V(/5*(VP)) > 0. If dtW — Qw < 0, then max w = 

max If;. 

Ft 

(2) Similarly, suppose that w' S H C(fdT) and div(V(/5*(VP)) < 0. If dtw' — Qw' > 0 , 

then min w' = minw'. 

Ut 

Proof. (1) This statement follows directly from Proposition 14.II bv taking w' = 0. 

(2) Similarly, this statement follows directly from Proposition 14.II bv taking w = 0. 

□ 


Now we are ready to prove the main results for the original operator L. 

Proofs of Theorem \1.4\ and Theorem \1.5[ 

Again, by the positivity and monotonicity of the exponential transformation u = exp(w), Theorem 
11.41 and Theorem 11.51 respectively follow from Proposition l4.1l and Proposition 14.21 □ 


5 Conclusion and outlook 

In this paper, we derive a generalised relativistic heat equation with external force helds and 
study comparison and maximum principles for this equation. The latter has been carried out 
using an logarithmic transformation which hinges on the gradient flow structure of the Boltzmann 
entropy with respect to an optimal transport associated to the general relativistic cost function. 
The derivation opens some interesting problems. The first possibility is to derive Eq. © use only 
optimal transport theory. It is not clear which energy functional should be used. The second open 
question is to rigorously prove that when c goes to infinity, Eq. ([3]) degenerates to the classical linear 
drift-diffusion equation. Research in this direction has recently been received lot of attention. 
Eor example, in mm. the author has proved that solutions of the relativistic heat equation 
without external force field indeed converge to that of the classical heat equation. We believe 
that the methods in these papers can be extended to include the external force fields. The third 
open problem is to derive Eq. ©, and in general flux-limited diffusions, from stochastic models. 
Connections between partial differential equations and stochastic processes often provide insightful 
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information about the geometrical structure and also offer new techniques for coarse-graining of the 
former, see e.g. [15111411131112] for more information in this direction. As mentioned in Remark 
13.31 the principle part of the operator Q is exactly the well-known mean curvature operator. 
Derivation of the mean-curvature flow using stochastic approach has been studied extensively in 
the literature, for instance in [23 [24]. It would be interesting to investigate whether one can 
derive Eq. and in general flux-limited diffusions, using such approaches. We leave these open 
problems for further research in the future. 
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